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Abstract 

The issue of developing effective and robust schemes to implement 
general hyperelastic constitutive models is addressed. To this end, special 
purpose functions are nsed to symbolically derive, evaluate, and automat- 
ically generate the associated FORTRAN code for the explicit forms of 
the corresponding stress function and material tangent stiffness tensors. 
These explicit forms are valid for the entire deformation range (i.e.,with 
both distinct and repeated principal-stretch values). The analytical form 
of these explicit expressions is given here for the case in which the strain- 
energy potential W is taken as a nonseparable polynomial function of the 
principle stretches. 


1 Introduction 

Recently, constitutive models of rubber hyperelasticity, using alternative rep- 
resentations in terms of the principal stretches as opposed to deformation in- 
variants, have become increasingly popular in nonlinear finite element analyses 
[1-3]. Two of our recent publications have discussed in detail the theoretical 
development [4] and symbolic and numeric implementation [5] of explicit forms 
for the second Piola Kirchhoff stress tensor and the material tangent stiffness 
tensor. These forms correspond to a class of Ogden type hyperelastic constitu- 
tive models, based on principal-stretch values. These explicit expressions are 
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significant since they are valid for the entire deformation range, even though 
the main constituents of the deformation tensor (i.e., principle values and asso- 
ciated eigenvectors) are, in general, neither uniquely defined nor continuously 
differentiable over the entire range. The two specific forms of the Ogden-type 
strain energy functions addressed in reference 4 encompass many of the popular 
representations currently in use for rubber materials. However, those functions 
were restricted to special f orms of nonseparable representations of the strain en- 
ergy density functions, with the restricted nonseparable form given in reference 
4, section 5, dealing with the important and practical case of incompressible 
and slightly compressible solids. To date, comparable treatments for the gen- 
era/ nonseparable forms of the models are not available in the literature. Indeed, 
it is the extension of our earlier results [4] and recent developments [5] to deal 
with this latter case that constitutes our main objective in the present paper. 

By cleverly applying symbolic manipulation packages so as to control ex- 
pression growth new constitutive theories can be developed and applied (e.g., 
finite element; see, [6] and [7]). Symbolic computation uses numbers, formulas, 
vectors, matrices, equations and the like to compute exact solutions; whereas 
numerical computation uses floating-point numbers to compute approximate 
solutions to problems of practical interest. Here, we will utilize three recently 
developed [5] special purpose symbolic functions (SDIFF, SDIFFEV, and TEM- 
PLATE) running under DOE MACSYMA [8]. These special purpose functions 
allow the derivation and automatic FORTRAN code generation of alternative 
generalized potential based constitutive models composed of principal values 
and their associated eigenvectors. 

This paper begins by reviewing highlights of our previous work in developing 
the theory of explicit forms [4] and implementing them symbolically and numer- 
ically [5]. Following this review, the results of the derivation of the generalized 
expressions for the second Piola Kirchhoff stress tensor Sij and the material 
moduli tensor Dijki are given. The paper then concludes with a discussion of 
the template files required to automatically generate the associated FORTRAN 
source code. 

2 Background 

The theoretical development of a singularity- free representation of principal 
value-based constitutive models has been discussed at length in reference [4], 
Here, we will confine our discussion to hyperelastic isotropic materials whose 
strain energy function W is taken to be a general function of the principal 
stretches, that is, 


W = W(X lt \ 3 ,Xa) (1) 

where A, , A2, A3 are the principal values of the right Cauchy-Green deformation 
tensor C,,. Denoting n, (*' = 1, 2, 3) to be the associated eigenvectors of C tJ , we 
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can define 


< 2 ) 

/=i 

where which is often referred to as the (orthogonal) eigenprojection oper- 
ator related to the associated eigenvectors of Cij is defined as 

N$P = njnj (3) 

Equation (2) is valid when all three eigenvalues (A,) are distinct. However, when 
two eigenvalues are the same (i.e, double coalescence, Ai ^ A 2 = A 3 = A), we 
have 

Cij = (Ai - \)NV + X6 {j (4) 

And for the case of triple coalescence (Ai = A 2 = A 3 = A), we have 




/=1 


(5) 


Similarly, through suitable manipulation of equations (2) and (4), explicit 
expressions for N$p in terms of Cij can be obtained for the case of three distinct 
eigenvalues, 

- (A,-A,)(A,-A,) KC * - (6) 

and for the case of double coalescence, 

-**«>■ <7) 

In the preceding equations r, s, and t represent any cyclic permutation of (1, 2, or 
3). These definitions will prove very useful in obtaining the pertinent singularity- 
free directional derivatives of both the potential W and the stress function 
Sij = Si j (Cij). 

The explicit singularity- free expressions for the second Piola Kirchhoff stress 
tensor Sij (Cij) are defined as 

Sy=2^sSy(Ctf) (8) 

Those for the material moduli tensor Dijki(Cij) can then be obtained by apply- 

ing the directional derivative formula to Su, that is, 


Dijtt = 



d 2 w 

dCijdCki 


Dijki(Cij) 


(9) 
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As a result, the explicit expressions for the tensors Si> (Cij) and D{j kiiCij ) 
can be obtained directly for the following three cases: case I - when all three 
eigenvalues are distinct; case II - when a single singularity is present (Ai / 
A 2 = A3 = A, i.e., double coalescence); or case III - when a double singularity is 
present (Ai ^ A 2 = A 3 = A, i.e., triple coalescence). 

The derivation and implementation process for the formulations described 
was recently automated [5] by constructing three special purpose functions (SD- 
IFF, SDIFFEV, and TEMPLATE), written at the MACSYMA command level, 
that can respectively, 

(1) Derive the explicit expressions for the stress tensor Si: (eqs. (8)) and 
material tensor Dijki (eqs.(9)), given three, one or no distinct eigenvalues 

(2) Evaluate symbolically the expressions generated by SDIFF for a given 
strain-energy function W 

(3) Evaluate the expressions generated by SDIFF and automatically generate 
(using the built-in MACSYMA function gentran) the associated FOR- 
TRAN code needed to evaluate the expressions numerically for a given 
potential function, W 

These three special purpose functions contain a list of built-in MACSYMA 
instructions (factor, expand, ev, ratsubst, diff, limit and for- loops, to 
name a few) arranged in a specific algorithmic order. Thus each special purpose 
function can be thought of as a macro command. 

3 Symbolic Derivation 

Let us begin by assuming that W is a nonseparable function of Ai, A 2 , and A 3 . 
For example, 

v 

W = ^[zn(Ai + A 2 4* A3) a * + 3 / n (AiA 2 + A 2 A3 + AsAi)^* + 2n(AiA 2 A3) 7m ] 

n = l 

As a result, when the special purpose function SDIFF is invoked, the scalar 
derivative of W with respect to each eigenvalue will no longer be a function of 
that eigenvalue only, as discussed in reference 5, but will instead be a function 
of all three eigenvalues, that is, 

«(i)(Ai,A 2 > A 3 ) = 2 ^“- 

Furthermore, in deriving the mixed second derivatives ( must also 

be taken into account in the procedure. To derive the generalized explicit ex- 
pression for the three cases, one need only issue the command SDIFF upon 
invoking MACSYMA, as shown here: 
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• Case I -three distinct eigenvalues (Ai ^ A2 ^ A3) 


SDIFF( 1 ) 

• Case II - double coalescence (Ai ^ A2 = A3 = A) 

SDIFF( 2 ) 

• Case III - triple coalescence (Ai = A2 = A3 = A) 

SDIFF( 3 ) 

Note that the resulting derived expressions have been manipulated and con- 
densed so as to streamline their reporting and to facilitate their comparison 


with previous work [ 4 ]. 

3.1 Results for Case I 

The explicit expression for the second Piola Kirchhoff stress tensor is 

Sij — ClCikCjcj + bC%j + c6ij ( 1 ®) 

where Sij is the second order identity tensor and a,b, and c are defined as 

a = — m[si (A2 — A3) + $2(^3 — ^1) + ®a(^i ” ^2)] (H) 

b = m[t,( Ai - A|) + s 2 (Xl - A l) + s 3 ( - A 2)] ( 12 ) 

c = — m[siA2A 3 (A2 — A 3 ) + 52A 3 Ai(A 3 — Ai) + $3^i^2(^i — A2)] ( 13 ) 


and where 

m = (Aj - A 2 )(A 2 - A 3 )(A 3 - AO (H) 

The explicit expression for the material moduli tensor Dijti(Cij) is 

Dijkt ~ aiP(Chj Cfj) + a 2 [P(Ch , Cij) + P(Cki y Cfj)] 

+ a 3 [Q(Cl6ij) + P{S hU Cfj)) + a 4 P(C kh Cij) ( 15 ) 

+ a 5 \Q{Cku j) + Q{S kh Cij )] 4 * 2 a$Iijki 

where two second order symmetric tensors P and Q have been introduced and 
are defined as 


P ljJb/ (G, H) = GikHjj + GuHjk ( 16 ) 
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( 17 ) 


(18) 


(19) 


Qijki(G, H) = GikHji + GijHjt + GjiHik + G jh Hu 

and the notation 

Iijkl = gftfcfy + tiutijk] 

Cfj = CtmCmj 

has been used in equation (15). Here the coefficients fli, 02 > are defined as 

3 3 3 

( 20 ) 
( 21 ) 
( 22 ) 


<*2 


a i=£’?»- + ][2 £ 

r=l r— 1 5=1, r^5 

3 3 3 

= £(Ar - h)vr - 2 £ £ (2Ji - A r - A,)£r 


r=l 

3 


’ r=l 5=l f r*5 
3 3 


a 3 = £x L + |£ £ (Ar + A.XA-Ar-A.tfr 

r=l r r=l 5 = l,r?fa 


a 4 = ][>- A r )V + |E £ ai-M(A-A,)f 

r=l ' ' 


3 3 


(23) 


' r=l 5=1, r^5 

3 3 


a 5 = X> + /3,?r( A~ Jl) )-l£ £ a2 + A r A,)(/i-A r -A,X„ (24) 

r=l r r= 1 j = l, r^J 

3 3 


r=l 


a 6 = ^(^) 2 77 r + (A r - /i)Pr + ^ h(h — A r — A 3 )^r, (25) 

A " r=l 5 = 1, r^5 


where 


and 


A*r = 


(Ar - A,)(A r - A t ) 


[s rr + (A, — A r )(^r + /*t) + ( At ~ A r )(/*r + Ps)] 
(A r — A,) 2 (A r — At ) 2 


£rs — 


S rs 


(A r - A a ) 2 (A r - A t )(A, - A t ) 

Ji — Ai + A2 + A3 
J 2 = A 1 A 2 4* A 2 A 3 + A 1 A 3 
h = A1A2A3 

Note that in the preceding expression the following differentiation notation has 
been introduced: 


6 


and 


s/ = si (Ai i A2, A3) — 2 


(26) 


Sjj — Sjj (Ai , A2, A3) — 


a«(Ai.A 2 ,A3) _ d 2 W 


d\j 


dXiXj 


(27) 


For example, « n = 2^r;s 2 2 = 2 ^t 5 s 33 = 2 aif; s i2 = s 2i - ^^’ 5l3 - 

s 3 i = ; s 2 3 = S32 = T-- A comparison of the preceding expressions 

and those obtained earlier for the two special Ogden-type strain energy lorms 
[4], shows that the expressions are identical except for the additional double- 
summation terms (containing the cross derivative terms) in the coefficients 
ai,a 2j ...a6 (see equations ( 20 )-( 25 )) comprising the material moduli tensor 
Dijkb Thus the previous work is now merely a special case of the present 
generalized expressions . 


3.2 Results for Case II 

In this case, a single singularity (Ai ^ A 2 = A3 = A) is analytically removed, 
thereby yielding 


with 


Sij = aCij + b&tj 


- S 2 

a = 


(Ai-A) 

t [$iA - s 2 Ai] 

(Ai - A) 


( 28 ) 

( 29 ) 

(30) 


and a reduced material moduli tensor 

Dijki = h P(C kl , Ch ) + h[Q(C kl , 6 {j ) + Q(S k i, Cij)] + b 3 I ijkl ( 31 ) 


where 


h = 


1 


(Ai - A 2 ) 3 


{(Ai — A 2 )[sh + s 2 2 — 2$i 2 ] — 2[si — s 2 ]} 


6 2 = ^ > ^ {(^2 — Ai)[A 2 su + Ais 22 - 2(Ai 4- A 2 )si 2 ] + (Ai 4- A 2 )[si - $ 2 ]} 

(Ax - * 2 ) 

&3 = 77 r~^r{(Al ~~ ^2)[^1 5 22 + AjSn — 2 A 1 A 2512 ] — 2AiA 2 [$i ~~ -S 2 ]} 

(Ai — A^ 

Again, in comparing the coefficients a and b, and, 61, i> 2) and 63 to those obtained 
in previous work [4], the only difference seen is the appearance of the cross 
derivative term ($i 2 ) in coefficients 61, & 2 , and 63. 
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3.3 Results for Case III 

Finally, in the case of a double singularity (A x = A 2 = A3 = A), the explicit 
expression for the stress tensor becomes 

Sij = s\(\)6ij (32) 

whereas the material moduli tensor becomes 

Dijkj = ( 33 ) 

These are identical to the previous results, as one would expect. 

The value of automating the foregoing derivation procedure is apparent in 
that not only does this special purpose function SDIFF relieve the user of the 
tedious manual derivation process, but it also ensures analytical accuracy. This 
was illustrated prior to the publication of reference 4, in that a number of errors 
in the hand derivation were detected, verified, and corrected. Also, because the 
derivation process needs to be executed only once, except for the evaluation of 
the scalar derivatives in equations (26) and (27) for each new definition of W, a 
second special purpose function, SDIFFEV, as described in [5], was developed. 
This function is used to symbolically evaluate the foregoing expressions. 

3.4 FORTRAN Code Generation 

The function TEMPLATE is similar to the function SDIFFEV in that both func- 
tions will evaluate the explicit expressions obtained from SDIFF. TEMPLATE, 
however, will automatically generate the associated FORTRAN source code 
needed to numerically evaluate the expressions for a given potential function 
W. Code generation is accomplished by utilizing the MACSYMA built-in func- 
tion gentran and a number of template files. The template files can be thought 
of as a framework for the FORTRAN generation of four subroutines (the main 
driving routine COMPSD and three subroutines, one each for case I , case II, 
and case III) and numerous functions. The template file for the main driving 
routine COMPSD is shown in appendix A. This subroutine is constructed for 
easy implementation into a finite element code. The input requirements are the 
strain tensor C m (denoted as emu) and its associated eigenvalues Ai ( A 2 , and A 3 
(denoted by gll, gl2, and gl3 respectively). The outputs are the stress tensor 
S n (denoted as s) and the material moduli tensor D nm (denoted as d). Here, 
n and m run from 1 to 6. Clearly, the only code generation required is that 
of subroutines COMPSD1, COMPSD2, and COMPSD3. Code generation is 
initiated by issuing the command gentranin, preceded by and followed by less 
than and greater than symbols, respectively. 

The subroutines COMPSD1, COMPSD2, and COMPSD3 are associated 
with case I (Ai ^ A 2 ^ A3), case II (Ai,A = A 2 = A 3 ), and case III (A = 
A x = A 2 = A 3 ), described in section 2.0. The template files corresponding to 
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these three cases are shown in appendices B,C and D, respectively. Note that in 
these routines, most of the FORTRAN code is automatically generated, since 
it pertains to the definition of coefficients a,b,c ; ai,a 2 , and the first and 

second scalar derivatives of the strain energy function W, (i.e., $i , $ 2 > $ 3 , $n > «22 > 
and S33). Also, the gentran commands are again preceded and followed by 
double inequality signs (that is, <C >>). All functions that are associated with a 
given case have been included in the corresponding appendix. As a result, with 
the appropriate template files, the FORTRAN source code associated with any 
general nonseparable or separable strain-energy function can easily be generated. 

4 Summary of Results 

Taken separately, the main constituents of the deformation tensor (i.e., the prin- 
cipal values and associated eigenvectors) are, in general, not uniquely defined 
and continuously differentiable functions. Careful consideration is thus called 
for in implementing constitutive models formulated in terms of these principal- 
strain measures. This difficulty was entirely bypassed by resorting to explicit 
symbolic derivations of appropriate forms of the material tangent-stiffness mar 
trices which are valid for the entire deformation range. Furthermore, to enhance 
effective utilization and implementation of the present results, automatic FOR- 
TRAN code generation of the present generalized explicit expressions was pur- 
sued and achieved. As a result, nonseparable forms dealing with the important 
practical case of incompressible and slightly compressible solids can easily be 
generated. Finally, the generic analytical forms of these explicit expressions have 
been given for three cases: (1) distinct eigenvalues, (2) one distinct eigenvalue, 
and (3) no distinct eigenvalues. 

In the future we will broaden our scope of application to include not only 
deformation constitutive models but also damage representations as well. An 
example that immediately comes to mind, where the above singularity-free rep- 
resentations will be important, is a maximum principle stress (or strain) damage 
formulation. Using this work as a building block, we can then envision moving 
to even more sophisticated damage formulations involving even higher tensorial 
representations. 
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APPENDIX A: Template File Associated With COMPSD 
The Main Driver Routine 


c ***************************************************** 
c This is the template subroutine to calculate 

c tensor S and D. inputs are eigenvalues gll,gl2,gl3, 

c and cmu(6) . emu is assumed to be engineering strain(e) , 

c e.g. the Cauchy-green deformation tensor cm(3,3) is related 

c to emu (6) in the following fashion: 

c cm(l , l)=cmu(l) , cm(2,2) = emu (2) , cm(3 ,3) =cmu(3) , 

c crau(4)=2*cm(l ,2) , cm(5)=2*cm(2,3) ,cmu(6)=2*cm(l ,3) . 

c The outputs are the second order tensor S(6) 

c and forth order tensor D(6,6) are related in the 

c following way: 

c S=D*C 

c S(l,l) = S(l) 

c S(2 ,2) = S(2) 

c S(3,3) * S(3) 

c S(l,2) * S(4) 

c S(2,3) = S(5) 

c S(3,l) = S(6) 

c C(l,l) = C(l) 

c C(2,2) = C(2) 

c C(3,3) - C(3) 

c C(1 ,2) = C(4) 

c C(2 ,3) = C(5) 

c C(3,l) = C(6) 

c 

subroutine compsd(gll ,gl2 ,gl3 , emu, s ,d) 
real*8 gll,gl2,gl3,ts(3,3) ,td(3,3,3,3) 
real*8 delt(3,3) ,delt4(3,3,3,3) ,s(6) ,d(6,6) 
real*8 cmu(6) ,cm(3,3) 
c 

c converts cmu(6) to matrix cm(3,3) in a way that 

c cm(l , 2)=cm(2 , 1) =cmu(4) , cm(2 ,3)=cm(3 ,2)=cmu(5) , 

c cm(l ,3)=cm(3, l)=cmu(6) . 

c 
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do 5 i=l,3 
do 5 j=l,3 

if (i.eq.j) then 
iq=i 

cm(i , j)=cmu(iq) 
else if (i.ne.j) then 
if ((i+j) .eq.3) iq=4 
if ( (i+j ) . eq . 4) iq=6 
if ((i+j) .eq. 5) iq=5 
cm ( i , j ) = emu ( i q ) / 2 
end if 
cont inue 

5 continue 
c 

c Initiates the second identity tensor delt(3,3) which 

c is a 2X2 identity matrix, 

c 

do 6 i=l,3 

delt(i,i)=1.0 

6 continue 
c 

c Computes the forth order identity tensor delt4(3,3,3,3) 

c by definition, 

c 

do 7 i=l ,3 
do 7 j=l ,3 

delt4(i, j ,i, j)=delt(i,i)*delt(j , j)+delt (i, j)*delt (j ,i) 
delt4(i , j , j , i)=delt4(i , j ,i,j) 

7 cont inue 
c 

c ******************************************************************* 
c For different eigenvalues gll,gl2,gl3 the computation 

c is different, easel is gll#gl2#gl3 call subroutine comsdl . 

c case2 is gl3=gl2#gll or gll=gl3#gl2 or gll=gl2#gl3 then 

c call subroutine compsd2. case3 is gll=gl2=gl3 call subroutine 

c compsd3 . 

c ******* ********************************************* **************** 
if ((gll .ne.gl2) .and. (gl2.ne.gl3) .and. (gll .ne.gl3)) then 
call compsdl(gll ,gl2,gl3,delt ,delt4,cm,ts ,td) 
else if ((gl2.eq.gl3) .and. (gll .ne.gl3)) then 
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call compsd2(gll ,gl2,delt ,delt4,an,ts ,td) 
else if ((gll . eq.gl2) .and. (gl3.ne.gl2)) then 
gll=gl3 

call compsd2(gll ,gl2,delt ,delt4,cm,ts ,td) 
else if ( (gll . eq.gl3) .and. (gl2 .ne.gl3)) then 
gll=gl2 
g 12=gl3 

call compsd2 (gll ,gl2,delt ,delt4,cm,ts,td) 

else 

call compsd3(gll ,delt ,delt4,ts,td) 
end if 

Rewrite the tensor ts(i,j) td(i , j ,k,l)to S(i) and D(i,j) 
respectively by using the symetric property, 
converts ts(3,3) s(6) and td(3,3,3,3) to D(6,6) 

do 8 i=l,3 
do 8 j=i,3 

if (i.eq.j) iq=i 
if (i.eq.l.and.j .eq.2) iq=4 
if (i .eq.2. and. j . eq.3) iq=5 
if (i.eq.l.and.j .eq.3) iq=6 
s(iq) =ts(i, j) 
continue 
continue 
do 9 i=l,3 
do 9 j=i,3 

d(i,j)=td(i,i, j , j) 

continue 
continue 
do 10 i-1,3 

d(i,4)=td(i,i,l,2)+td(i > i,2,l) 

d(i,5)=td(i,i,2,3)+td(i,i,3,2) 

d(i,6)=td(i,i,3,l)+td(i > i,l,3) 

continue 

d(4,4)=(td(l ,2, 1 ,2)+td(l ,2 ,2 ,l)+td(2,l ,1 ,2)+td(2 , 1 ,2,l))/2 
d(4, 5)=(td(l ,2,2 ,3)+td(l ,2 ,3,2)+td(2,l ,2,3)+td(2 ,1 ,3,2))/2 
d(4,6)=(td(l ,2, 1 ,3)+td(l ,2 ,3 , l)+td(2, 1 ,1 ,3)+td(2 , 1 ,3 , 1) ) /2 
d(5,5)=(td(2,3,2,3)+td(2,3,3,2)+td(3,2,2,3)+td(3,2,3,2))/2 
d(5,6)=(td(2,3,l,3)+td(2,3,3,l)+td(3,2,l,3)+td(3,2,3,l))/2 



d(6,6)=(td(3,l ,1 ,3)+td(3,l ,3, l)+td(l,3,l ,3)+td(l ,3,3,l))/2 . 
do 11 i = 1,6 
do 11 j = 1,6 

d(i , j) = d(j ,i) 

11 continue 

c 

c prints out the inputs gll ,gl2 ,gl3,cmu(6) and outputs S and D 

c 

print*, ’gll*', gll 

print*, ’gl2=’, gl2 

print*, ’gl3=’, gl3 

print*, ’Input tensor C(6) : ’ 

print*, (cmu(i) , i = 1,6) 

print* , "second order tensor S(6):" 

print*, (s(i), i*l,6) 

print*, "The forth order tensor D(6,6):" 
do 101 i=l ,6 

print*, (d(i, j) , j=l,6) 

101 continue 

return 
end 
c 

subrout ine comp sdl (gl 1 , gl2 , gl3 , delt , delt4 , cm , t s , t d) 

« 

gentranin("casel .tem")$ 

» 

subroutine compsd2 (gll ,gl2 , delt , delt4 , cm, t s , td) 

« 

gentranin("case2 .tem")$ 

» 

subroutine compsd3(gll ,delt ,delt4,ts,td) 

« 

gentranin("case3 .tem")$ 

» 
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c 

c This subroutine computes P and Q forth order tensors 

c which we define in tensor D. 

c 

subroutine pqcom(cml,cm2,p,q) 
real*8 cml(3,3) ,cm2(3,3) , p(3,3,3,3) ,q(3,3,3,3) 
do 100 i=l ,3 
do 100 j-1,3 
do 100 k=l ,3 
do 100 1=1,3 

p(i,j ,k,l)=cml(i,k)*cm2(j ,l)+cml(i,l)*cm2(j ,k) 
q(i , j ,k,l)=p(i, j ,k,l)+cml(j ,l)*cm2(i,k)+cml(j ,k)*cm2(i,l) 
100 continue 

return 
end 
c 

c This subroutine computes matrix product cmXcm. 

c 

subroutine product (mat l,cmm) 
real*8 mat 1(3, 3) ,cmm(3,3) ,sum 
do 25 i-1,3 
do 25 j=l ,3 
s\im=0 . 0 
do 26 k=l ,3 

s\un=sxnn+mat 1 (i , k) *mat 1 (k , j ) 

26 continue 

cmm(i , j)=sum 
25 continue 
return 
end 
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APPENDIX B: Template File Associated With C0MPSD1 
Valid For Three Distinct Eigenvalues 


real*8 gll ,gl2,gl3,ts(3,3) ,td(3,3,3,3) 
real*8 an(3,3) ,delt(3,3) ,delt4(3,3,3,3) ,p(3,3,3,3) 
real*8 q(3,3,3,3) ,cmm(3,3) , pi (3, 3, 3, 3) ,p2l(3,3,3,3) 
real*8 p3l(3,3,3,3) ,qll(3,3,3,3) ,ql2(3,3,3,3) ,p22(3,3,3,3) 
real*8 q2l(3,3,3,3) ,q22(3,3,3,3) ,a,b,c,al,a2,a3,a4,a5,a6 
c 

c Obtains cmm(3 ,3)=cm(3 ,3)*cm(3,3) from subroutine product 

c 

call product (cm, cmm) 

c Uses the formula we derived in code to compute second order 

c tensor ts(3,3). 

« 

gent ran (for i : 1 thru 3 do 
(for j : 1 thru 3 do 

(ts [i , j] :a(gll ,gl2 ,gl3) *cmm[i , j] +b (gll ,gl2 ,gl3) 

*cm[i, j]+c(gll,gl2,gl3)*delt[i,j])))$ 

» 

c Call subroutine to compute all the functions we defined 

c when we derived forth order tenosor td, namely P(i,j,k,l) 

c and Q(i,j,k,l) which are the functions of cm(3,3) and 

c the matrix product cmm(3,3). 

c 

call pqcom(cmm,cmn > pl ,q) 
call pqcom(cmm,cm,p21 ,q) 
call pqcom(cm,cmm,p22 ,q) 
call pqcom(cm,cm,p31 ,q) 
call pqcom(cmm,delt ,p,qll) 
call pqcom(delt ,cmm,p,ql2) 
call pqcom(cm,delt ,p,q21) 
call pqcom(delt ,cm,p,q22) 
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c 

c Computes forth order tensor td(i,j,k,l) 

c 

« 

gentran(for i : 1 thru 3 do 
(for j : 1 thru 3 do 
(for k:l thru 3 do 
(for 1:1 thru 3 do 

(td[i, j ,k,l] :al(gll ,gl2,gl3)*pl[i, j ,k,l]+a2(gll ,gl2,gl3) 
*(p21[i,j ,k,l]+p22[i , j , k , 1] ) +a4 (gll , gl2 , gl3) *p3l[i,j ,k,l] 
+a3 (gll , gl2 , gl3) * (ql 1 [i , j , k , 1] +ql2 [i , j , k , 1] ) 

+a5 (gll ,gl2 , gl3) * (q21 [i , j ,k,l]+q22[i,j ,k,l]) 

+a6 (gll ,gl2 , gl3) *delt4 [i , j ,k , 1] ) ) ) ) ) $ 

» 

return 

end 

c 

c a,b,c,al-a6 are the coefficients we derived in code, 

c 

« 

gentran(a(gll ,gl2,gl3) :=block (type (function, a) , 

type("real*8" ,gll ,gl2 ,gl3) , 
type ("real *8" ,a, si , s2 ,s3) , 
a : eval (ta) ) ) $ 

» 

« 

gentran(b(gll,gl2,gl3) :=block (type (funct ion, b) , 

type C'real*8" ,b , gll ,gl2 , gl3) , 
type("read*8" ,sl ,s2 ,s3) , 
b : eval (tb) ) ) $ 

» 

« 

gentran(c(gll ,gl2,gl3) : =block(type (funct ion, c) , 

type ("real *8" ,c,gll ,gl2 ,gl3) , 
type("real*8" ,sl ,s2 ,s3) , 
c:eval(tc)))$ 

» 
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« 

gentran(al (gll ,gl2 ,gl3) :=block(type (function, al) , 

type("real*8" ,al ,gll ,gl2 ,gl3) , 
type ( "real*8 n , si , s2 , s3 , s 11 , s22 , s33 , s21 , s32 , s31 ) , 
al :eval(al)))$ 

» 

« 

gentran(a2(gll,gl2,gl3) :=block(type(f unction, a2) , 
type ("real*8",a2, gll, gl2,gl3) , 

type(' , real*8" , si , s2 , s3 , sii , s22 , s33 , s21 , s32 , s31) , 
a2:eval(a2)))$ 

» 

« 

gentran(a3(gll,gl2,gl3) :=block (type (f unction, a3) , 
type("real*8",a3,gll,gl2,gl3) , 

type ( "real*8" , si , s2 , s3 , sll , s22 , s33 , s21 , s32 , s31) , 
a3 : eval (a3) ) ) $ 

» 

« 


gent ran (a4 (gll ,gl2,gl3) : =block (type (function, a4) , 
type("real*8",a4,gll,gl2,gl3) , 

type("real*8" , si ,s2,s3,sll , s22 ,s33, s21 ,s32 ,s31) , 
a4:eval(a4)))$ 

» 

« 

gentran(a5(gll,gl2,gl3) :=block (type (funct ion, a5) , 
type("real*8" ,a5,gll,gl2,gl3) , 

type("real*8" ,sl ,s2,s3,sll ,s22,s33, s21 ,s32, s31) , 
a5:eval(a5)))$ 

» 

« 

gentran(a6(gll,gl2,gl3) : =block (type (f unction, a6) , 
type("real*8" ,a6,gll,gl2,gl3) , 

type("real*8" , si ,s2,s3,sll ,s22 ,s33, s21 ,s32 , s31) , 
a6:eval(a6)))$ 


18 



c 

c The sl,s2,s3,sll,s22,s33,s21,s32,s31 axe derivatives of V 

c 

function slCgll ,gl2,gl3) 

<<cut(var) ;» 

« 

gentran(type( "real*8" ,si ,gll ,gl2 ,gl3) , 

si :2*eval(diff (w, ’gll,l)))$ 

» 

return 

end 

c 

f unct ion s2 (gll , gl2 , gl3) 

<<cut(var) ;» 

« 

gentran(type( "real*8",s2,gll,gl2,gl3) , 

s2 : 2*eval(diff (w, ’gl2,l)))$ 

» 

return 

end 

c 

function s3(gll ,gl2,gl3) 

<<cut (var) ; » 

« 

gentran(type( "real*8" , s3,gll ,gl2,gl3) , 

s3:2*eval(diff (w, ^13,1)))$ 

» 

return 

end 

c 

function sll (gll ,gl2 ,gl3) 

<<cut (var ) ; » 

« 

gentran(type( "real*8" ,sll ,gll,gl2,gl3) , 

sll : 2*eval(diff (w, ’gll, 2)))$ 

» 

return 

end 
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function s22(gll ,gl2,gl3) 

<<cut (var) ; » 

« 

gentran(type("real*8" ,s22,gll ,gl2,gl3) , 

s22:2*eval(diff (w, ’gl2,2)))$ 

» 

return 

end 

c 

function s33(gll ,gl2,gl3) 

<<cut (var) ;>> 

« 

gent ran (type ( "real*8" , s33,gll,gl2,gl3) , 

s33 : 2*eval (dif f (*,^13,2)))$ 


» 

return 

end 

c 

function s2l(gll ,gl2,gl3) 

<<cut (var) ;» 

« 

gentran(type( "real*8" ,s2i ,gll ,gl2,gl3) , 

s21 : 2*eval (dif f (w , J gl2 , 1 , * gll , 1) ) ) $ 

» 

return 

end 

c 

function s31(gll,gl2,gl3) 

<<cut(var) ;» 

« 

gentran(type( "real*8" ,s31 ,gll ,gl2,gl3) , 

s31:2*eval(diff(w, , gl3, l,*gll,l)))t 

» 

return 

end 
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c 


function s32(gll ,gl2,gl3) 
«cut (var) ; » 


« 


» 


gentran(type( "real *8" , s32,gll,gl2,gl3) , 

s32 : 2*eval (dif f (w , » gl3 , 1 , » gl2 , 1 ) ) ) $ 


return 

end 


21 



APPENDIX C: Template File Associated With C0MPSD2 
Valid For Double Coalesence Case 

real*8 gll ,gl2,ts(3,3) ,td(3,3,3,3) 
real*8 cm(3,3) ,delt(3,3) ,delt4(3,3,3,3) , pi (3, 3,3, 3) 
real*8 q2(3 ,3,3 ,3) ,ql (3, 3, 3, 3) ,p(3,3,3,3) ,q(3,3,3,3) 
real*8 bl ,b2,b3, abar.bbar 
c 

c Computes second order tensor ts(i,j) based on the formula 

c derived in code, 

c 

« 

gentran(for i:l thru 3 do 
(for j : 1 thru 3 do 

(t s [i , j] : abar (gll ,gl2) *cm[i , j] +bbar (gll ,gl2) *delt [i , j] ) ) ) $ 

» 

c 

c Call subroutine to get P, Q which are defined in code, 

c 

call pqcom(cm,cm,pl ,q) 
call pqcom(cm,delt ,p,ql) 
call pqcom(delt ,cm,p,q2) 
c 

c Computes tensor td(i,j,k,l). 

c 

« 

gent ran (for i:l thru 3 do 
(for j : 1 thru 3 do 
(for k:l thru 3 do 
(for 1 : 1 thru 3 do 

(td [i , j , k , 1] :bl (gll , gl2) *pl [i , j ,k , 1] +b2 (gll , gl2) * 

(ql [i , j ,k,l] +q2 [i , j ,k, 1] )+b3(gll ,gl2)*delt4[i , j ,k,l] )))))$ 

» 

return 

end 
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c 

c abar,bbar are bl, b2, b3 functions derived in code, 

c 

« 

gent ran (abarCgll ,gl2) :=block( type (funct ion, abar) , 

type("real*8",abar,gll,gl2) , 
type("real*8" , ssl,ss2), 
abar : eval (abar) ) ) $ 

» 

« 

gentran(bbar(gll ,gl2) :=b lock (type (funct ion, bbar) , 

type("real*8" ,bbar,gll ,gl2) , 
type("real*8" , ssl,ss2), 
bbar : eval (bbar) ) ) $ 

» 

« 

gentran(bl(gll ,gl2) :=block(type(function,bl) , 

type("real*8",bl,gll,gl2) , 

type ("real *8" , ssl , ss2 , ssll , ss22 ,ss21) , 

bl :eval(bl)))$ 

» 

« 

gentran(b2(gll,gl2) :=block (type (funct ion, b2) , 

type("real*8" ,b2,gll,gl2) , 
type("real*8" , ssl , ss2 , ssll , ss22 ,ss21) , 
b2:eval(b2)))$ 

» 

« 

gentran(b3(gll,gl2) :=block(type(function,b3) , 

type("real*8",b3,gll,gl2) , 
type("real*8" , ssl , ss2 , ssll , ss22 ,ss21) , 
b3:eval (b3)))$ 

» 

« 

new: sub st( [^13=^12] ,w)$ 

» 
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c 

c ssl,ss2,ssll,ss22,ss21 are derivatives of W. 

c 

function ssl(gll,gl2) 

<<cut (var) ;» 

« 

gentran (type ( "real*8" , ssl , gll , gl2) , 

ssl :2*eval(diff (new, ’gll ,1)))$ 

» 

return 

end 

c 

function ss2(gll,gl2) 

<<cut (var) ;» 

« 

gentran (type("real*8" ,ss2,gll ,gl2) , 

ss2:2*eval(diff (new, *gl2,l)))$ 

» 

return 

end 

c 

function ssll(gll,gl2) 

«cut (var) ; » 

« 

gentran(type("real*8" ,ssll,gll,gl2) , 

ssll :2*eval(diff (new, ’gll, 2)))$ 

» 

return 

end 

c 

function ss21 (gll ,gl2) 

<<cut (var) ; » 

« 

gentran (type ("real*8" ,ss21 ,gll ,gl2) , 

ss21 :2*eval(diff (new, ’gl2,l, ’gll,l)))$ 

» 

return 

end 
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c 

function ss22(gll ,gl2) 

«cut(var) ;» 

« 

gentran(type("real*8" , ss22 ,gll ,gl2) , 

ss22:2*eval(diff (neww, ’gl2,2)))$ 

» 

return 

end 
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APPENDIX D: Template File Associated With C0MPSD3 
Valid For The Triple Coalesence Case 


c 

real*8 gll,ts(3,3) ,td(3,3,3,3) ,delt(3,3) ,delt4(3,3,3,3) 
real*8 ccl, abbar 
c 

« 

gentran(for i:l thru 3 do 
(for j : 1 thru 3 do 
(ts[i, j] :abbar(gll)*delt [i, j])))$ 

» 

« 

gentran(for i:l thru 3 do 
(for j : 1 thru 3 do 
(for k:l thru 3 do 
(for 1 : 1 thru 3 do 

(td[i,j ,k,l] : cel (gll) *delt4[i , j ,k,l])))))$ 

» 

return 

end 

« 

gent ran (abbar (gll) : =b lock (type (function, abbar) , 

type("real*8" , abbar,gll), 
abbar : eval (abbar) ) ) $ 

» 

« 

gentran(ccl(gll) :=block(type(function,ccl) , 

type("real*8" , ccl,gll), 
ccl :eval(ccl)))$ 

» 

« 

m: subst ( [ , gl3 =, gll , ^12=^11] ,w)$ 

» 
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c 


« 


» 


c 


« 


» 


function sssl(gll) 

«cut (var) ; » 

gentran(type("real*8" ,sssl,gll) , 
sssl :2*eval(diff (www, ’gll ,1)))$ 

return 

end 

function sssll(gll) 

<<cut(var) ;» 

gentran(type("real*8" ,sssll ,gll) , 

sssll :2*eval(diff (www, , gll t 2)))$ 


return 

end 
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